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ABSTRACT 


In this report a method is proposed to treat the problem of 
steady, two-dimensional, laminar, incompressible, high Reynolds number 
separated flow past thin airfoils. An integral form of the boundary- 
layer equations with interaction is used and the interaction between the 
inviscid and viscous flowfields is provided for by use of a thin-airfoil 
integral. A detailed documentation of the attempts at obtaining a 
solution is presented. Even though these attempts were essentially 
unsuccessful, they help to suggest paths for future research on this 
difficult problem. In addition, a useful survey of the current state- 
of-the-art of problems involving viscous-inviscid interactions in flow- 
fields with separation is given. 
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CHAPTER I 


INTRODUCTION 

The problem under consideration is the theoretical study 
of laminar, steady, incompressible, high Reynolds number separated 
flow past two-dimensional airfoils. The quantitative description of 
flowfields containing separation is one of the major unsolved problems 
of fluid mechanics. To appreciate the difficulty involved in the 
proposed research, consider first the laminar, incompressible, high 
Reynolds number flow along a rigid surface in the absence of regions 
of separated flow. 

Prandtl demonstrated that the effects of viscosity can be 
considered to be confined to a narrow region - with thickness proportional 
to the inverse square root of the Reynolds number - close to the 
surface. He proposed the boundary-layer equations as an accurate 
approximation to the Navier-Stokes equations for flows with large 
values of the Reynolds number. To first order in the inverse square 
root of the Reynolds number, the pressure distribution does not vary 
across the boundary layer and may be obtained by a computation of the 
potential flow past the surface. The boundary-! ay er equations may 
then be integrated with this pressure distribution as a known quantity - 
various accurate numerical techniques are available, for example. 

Van Dyke 1 has shown how corrections may be added to the first-order 
solution to account for streamline displacement and surface curvature 
effects. 
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In the event that regions of separated flow are present 

along the surface of the body, the viscous layer thickens rapidly 

and first-order perturbations may be introduced into the outer 

potential flow field. This violates the result of classical boundary- 

layer theory that the outer flow to first order is independent of 

the boundary layer. With this breakdown of boundary-layer theory in 

the region near the separation point, it would seem that the appropriate 

path to follow to compute the flowfield with separation present is to 

integrate the complete Navier-Stokes equations. 

At the time this research was initiated, it was felt that 

a numerical solution of the Navier-Stokes equations for this problem 

was beyond the state-of-the-art and up to the time of the writing of 

this report, this is still the case. Navier-Stokes solutions to the 

high Reynolds number separated flow situation are still not available 

although much progress is being made in the development of suitable 

2 3 

numerical techniques. Recently, Davis and Davis and Werle have 

developed a numerical scheme for integrating the Navier-Stokes equations 

for flow past the parabola and paraboloid of revolution at all values 

of the Reynolds number. It remains to be shown whether this technique 

4 

will prove suitable for the separated flow problem. Jacob uses a 
boundary layer - inviscid outer flow iteration scheme to treat airfoil 
separated flow problems but uses a crude dead-water model for the 
separated flow and wake regions. The results compare well with available 
experiments. Similar treatment for the circular cylinder problem is 

5 

presented in Bluston and Paulson . 
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An appropriate technique for the solution of this problem 

will have to allow for interaction between the outer inviscid flow 

and the inner viscous layer. Strong interaction problems of this type 

where the viscous and inviscid layers develop simultaneously have been 

successfully treated for the case of a supersonic outer stream by Lester 

Lees and co-workers at the California Institute of Technology using 

a method known as viscous-invi>scid interaction theory. Some illustrative 

6 7 

examples of this method are presented in Klineberg and Lees and Alber . 

In the latter, Alber outlines an approach for the solution of the 
incompressible turbulent wake problem. 

The essential features of the viscous-inviscid interaction 
theory are as follows. It is assumed that the boundary-layer equations 
adequately describe the important features of the flow in the viscous 
region. Their use is modified, however, since the outer flow pressure 
gradient is not known in advance but is calculated simultaneously with 
the variables in the viscous region. The partial differential equations 
valid in the viscous layer are integrated across the layer and are thereby 
reduced to a set of ordinary differential equations. An equation is 
introduced to link the outer inviscid and inner viscous flowfields, 
and the complete set of equations is integrated with the inviscid 
surface speed as one of the unknowns. For the case of the supersonic 
outer stream, the linking equation usually is of the form of a Prandtl- 
Meyer relation. 

At the initiation of this research, it was proposed to extend 
the method of viscous-inviscid interactions to allow for the solution 
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of the incompressible flow problem under consideration. The main 
task was then to develop an appropriate equation or technique to 
link the inviscid and viscous flow regions. This task is complicated 
since the inviscid flowfield is governed by an elliptic partial 
differential equation (Laplace's equation) for the incompressible case 
so that the surface speed at any body location is dependent on the 
complete displacement thickness distribution. It was clear that to 
discover the appropriateness of this technique for the solution of 
this problem it would be useful to choose a flow configuration which 
would allow for simplifying approximations in the treatment of the 
inviscid flowfield. A thin symmetric Joukowski airfoil seemed an 
ideal choice since simplifications could be introduced to take advantage 
of both the symmetry and the small flow disturbances. A thin-airfoil 
integral could then be used as the linking equation and the need to 
compute the detailed inviscid flowfield is removed. 

In the following chapters, the method is developed and 
several attempts at solution are discussed. Since a solution was 
hot obtained, it was felt most appropriate to present the results 
of the research in an essentially chronological sequence to allow 
the reader to determine the approaches attempted and their subsequent 
outcomes . 

The problem for separated flow past biconvex airfoils has 
recently been solved for both subsonic and transonic flowfields by 

Q 

Klineberg and Steger of NASA-Ames Research Center. They found it 



necessary to compute the complete inviscid flowfield using the 
finite-difference relaxation technique developed by Steger and 
Klineberg and to make use of an interactive computer graphics 
system. Also, attempts are currently underway to solve the interaction 
problem without resorting to an integral representation for the 
viscous region. Klemp and Acrivos^ and others have developed 
numerical techniques for the integration of the boundary-layer 
equations through a reverse flow region and Erdos, Baronti and 
Elzweig^ discuss the inclusion of such a technique in a method for 
the solution of a transonic interaction problem. 

At present, research is underway to study interaction problems 
for flows ranging from incompressible to hypersonic and for both 
the laminar and turbulent cases. In this report, in addition to a . 
description of the present research (Chapters II and III), a detailed 
discussion is given of the results of Klineberg and Steger 8 
(Chapter IV) and a survey is presented of the current state-of-the- 
art in the area of viscous-inviscid interactions (Chapter V). 
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CHAPTER II 


MATHEMATICAL FORMULATION OF PROBLEM 
A. Viscous Region 

Consider the two-dimensional steady laminar flow of an 
incompressible fluid at high Reynolds number. Classical boundary- 
layer theory would require a consideration of the continuity and 
streamwise momentum equations. For the formulation of this interaction 
problem, it is necessary to include the mechanical energy equation 
which is obtained by multiplying the streamwise momentum equation by 
the streamwise velocity component. The coordinate system used is 
shown on Figure 1. s and n are distances measured along and normal 
to the surface, respectively, u and v are the velocity components in 
the s and n directions, p is the pressure, p is the density and y is 
the dynamic viscosity coefficient. 

The governing partial differential equations are 
Continuity: = 0 (1) 


Streamwise Momentum: 


0 ( u I” +v M) = .dE+ 
p ' u as v an' ds y an-' 


Mechanical Energy: 


/ o au . au\ dp . a 2 u 

c(u S? + uv 35-) ' ar + ^ -5H2 


The equations are now integrated across 
from the surface n = 0 to the boundary-layer edge 
u = v = 0. At the edge, u = u e (s), v = v g (s) and 
yau/3n = 0. The integrated equations are 


( 2 ) 

(3) 

the boundary layer, 
n = 6. At the surface, 
the shear stress 
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r 



( 4 ) 


tan e = 


d6 

ds 



dJtnu, 
e 

ds 


d 

cHF 


u|8 * 


* 

<5 U, 



y 


3 LI I 

an ■ n=0 


/p 


( 5 ) 



where v = y/p, the kinematic viscosity, e is the streamline direction 
at the boundary-layer edge, and 6 , 9 and 0 are the displacement, 
momentum and energy thicknesses given by 



It is convenient to introduce the following quantities which 
depend on the streamwise velocity profile. 

H = e/6* 


* * 

J = e /6 
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One profile shape parameter, say a(s), is chosen and the 
other profile quantities in Equation (8) will be determined as functions 
of it. Equations (4-6) can then be written: 


* * Hn 

d6 _ Z6_ au e 
3sF “ u 0 d~s 



= tan 0 

(9) 

H a? ♦^-<2 h+D 

du e 

. * dH da 

6 da" ds" 

= if-. 

U e 6* 

(10) 

, d6* . 3J6* du e 

J ar u 0 ds 

e 

* 

+ 6 

dJ dH da 
HIT da ds 

_ vR 

U 6* 

e 

(ID 


The unknowns in these three ordinary differential equations 

are 6*, u e> 0 and a(s). This set of equations is known as the strong 

interaction equations when u fi is considered unknown. The outer inviscid 

flow is then coupled to the inner viscous flow. 

The flowfield can be divided into four distinct regions in 

terms of the nature of the velocity profiles. They are illustrated in 

Figure 2. Region I, the body-attached flow region, extends from the 

leading edge to the separation point. Region II, the body-reverse 

flow region, extends from the separation point to the trailing edge. 

Region III, the wake-reverse flow region, extends from the trailing 

edge to the rear stagnation point of the wake. Region IV, the wake 

forward flow region, extends from the wake rear stagnation point to 

downstream infinity. Following Klineberg and Lees 6 and Klineberg, Kubota 
1 ? 

and Lees , a(s) is defined in the four regions as 
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3(U/Uj 


Region I 

a = 

a (n/6 ) n=0 

Region II 

a = 

V u/u e = 0 

Region III 

a = 

r— i 
V*=o 

Region IV 

a = 

r— l 

L ujn=0 


a(s) is therefore the wall shear stress, the thickness of the reverse 

flow region and the velocity ratios on the dividing streamline ip * 0 

and the wake center line, respectively. 

Similar solutions to the classical bound ary -layer equations 

are available for profiles characteristic of each of the regions. It 

is assumed that the functional relationships between the profile quantities 

in Equation (8) and a(s) are the same as those for the similar solutions. 

It is noted, however, that no dependence between profile shape and local 

pressure gradient is assumed. Extensive curve fitting has led to accurate 

polynomial representation of the profile quantities as functions of a(s). 

These are presented in Appendix A for each region. The results for 

the profiles in regions I and II are those in Klineberg and Lees 6 and 

for profiles in regions III and IV are those in Klineberg, Kubota and 

Lees . This choice of a(s) is made because of the availability of 

the above curve fits. A better choice might be H since it has one 

8 

definition in the flowfield r Klineberg and Steger use H in their analysis. 

* 

Boundary-layer theory predicts that tan e, 5 and 6 are of 
order v^. To work with variables of order one, the following scaling 
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is introduced. 

(tan ©) ' = v "^ 2 tan © 

(.*)■• =v-' /2 / 

«• -v' 1 ' 2 . 

The new variables are substituted into Equations (9-11) and the primes 
are dropped. To put the equations in a suitable form for numerical 
integration, Kramer's Rule is used to isolate the first derivatives. 
The equations become 

d 6 * = Dj ^e = D 2 da _ D 3 
<Js ” D * ds ~D’ds “ D 


(13) 


(14) 


where 


1 - tan 0 [(2H+1) - 3J] + Z (P gjj- - R)/u p 6 


3h 


D 2 " (P gjj- - R)/5* 2 - u e tan 0 (H^ - J)/6* 


D 3 = [R(2H+1+ZH) - PJ(3+Z)]/(u e 6 * 2 gj) + j tan 0 (H-l )/ 6 * jjjj- 


D = 3H - (2H+1+ZH) - J(3+Z) 


B. The Interaction Equation 

There are now three viscous flow equations and four unknowns, 
one of them being u g , the inviscid streamwise velocity component at the 
boundary-layer edge. It is necessary to introduce an equation to close 
the system and this equation is the Important one which provides for 
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the interaction between the viscous boundary layer and the outer inviscid 
flow. For problems with a supersonic outer stream, a Prandtl -Meyer 
relation is usually used to link 0 and u^. The problem becomes complicated 
in incompressible flow because the inviscid equations are elliptic and 
therefore global in nature. The local boundary-layer edge velocity 
depends on the displacement thickness distribution in the complete flow- 
field. 

Consider the flow schematic shown in Figure 3 where is the 

free stream speed, x and y are along and normal to the chordline. Let 

U and V be the velocity components in the x and y directions, respectively. 

With viscous effects included, the effective body shape is obtained by 

* 

adding the displacement thickness to the body. Note that 6 is measured 
normal to the body surface. Let y * f(x) represent the body and let 
its geometry be such that thin-airfoil theory approximations are valid. 

The effective body shape is represented by a distribution 
of sources along the x - axis of strength q(x) per unit length. The 
velocity potential for the irrotational flow is then 

4>(x,y) ’ x + 57 f q(§) log [(x-§) 2 + y 2 ] 1/2 d§ 

'0 

where 

q(x) = 2 U_ ^ [f(x) ♦ .>/* «*/(! ♦ (|i } 2 ) 1/2 1 

To the order of approximation of thin-airfoil theory, (df/dx) 2 «l and 
can be neglected in the above equation. Now, 


U 



and 


U = ^- = U 

9 X 



. 1/2 dS 

+ V T— 

dx 


(§)] ( x -§ V +? d§ 




, 1/2 d6 

+ v Hx 


(§)] 


W 


ds 


1/2 _ v(x,f+v 1/2 S) 

v (tan e) invisc , d - if ' 


° f f Eg|<§) -«■ v 1/2 jg- ]d S 

* * dX dX fx^l2+tf+u ,/Z Al2 


(x-§ ) 2 +(f+V o)2 

U 

Note that (tan 0)-j nvl * sc *id ** measured with respect to x while tan 0 is 
measured with respect to s or the local body slope. Therefore, 


1/2 


tan 0 = v ,/2 (tan e) 1nv1scid -■ g- 
and finally, the equation linking the flowfields is 


df 


tan 0 


f(x)+v 1/2 6(x) 

ir 



1/2 


$■*(«)] 


d§ 

(x-§)2+(f+v 1/Z S)2 


- v' 1 / 2 *$) (15) 

C. Matching Conditions 

For the viscous calculation, there are four distinct regions with 

different definitions of profile parameter a(s) and different polynomial 

relationships between the profile quantities. It is therefore necessary to 

enforce matching conditions as the solution progresses from one region into 

another. Physical considerations lead to the condition that displacement 
* 

thickness 5 and inviscid speed u 0 are continuous across region boundaries. 
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Since a(s) is defined differently in each region, conditions on it 

are considered separately. In region I, a(s) is the wall shear stress 

and the region ends at the separation point a = 0. In region II, a 

is the thickness of the reverse flow region and a = 0 at the start of 

the region. In region III, a is the dividing streamline velocity ratio. 

At the trailing edge, let a(II) be the value in region II and a(III) 

12 

be the value in region III. Klineberg, Kubota and Lees have developed 
a curve fit for a(III) as a function of a(II) and it is used as the 
trailing edge matching condition. 

a(III) = .05496 a(II) + .29702 a( II) 2 (16) 

+ 11.19943 a( II) 3 - 21.66525 a(II) 4 + 10.06854 a( II) 5 
At the end of region III, the wake rear stagnation point, a =0. At 
the start of region IV, a is the wake centerline velocity ratio and is also zero. 

D. Airfoil and Initial Potential Flow Solution 

An airfoil must be selected such that separation will occur 
but also such that thin-airfoil approximations will be appropriate. A 
symmetric airfoil at zero angle of attack is chosen so that the difficulties 
involved with correctly determining the circulation do not add to the 
complexity of the problem. For ease in calculation of the initial potential 
flow, a Joukowski airfoil section is used. 

13 

Van Dyke states that to second order a Joukowski airfoil 
is represented by 

y = e(2-x)(2x>-x 2 ) 1/2 (17) 

where the chord is 2 and the maximum thickness ratio is 3^ 2 e/4. For 
the 10% thickness ratio used in the analysis, e = .077. The surface speed 
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u g can easily be determined by use of conformal mapping. The 
Joukowski transformation maps the airfoil into a circular cylinder. 
To obtain u e (s), it is necessary to compute s(x). The airfoil, 
its slope and surface speed are given in Figure 4. 
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Chapter III 


METHODS OF SOLUTION 

The mathematical problem to be solved is given in Equations 

(14) and (15) with the appropriate polynomial representations of the 

profile quantities (Appendix A) and the matching conditions. The 

4 

Reynolds number considered is of the order of 10 . 

A. Karman-Pohl hausen on Airfoil Forward Portion 

Initially, it was felt that the upstream influence of the 
separated flow region would be of limited extent. Therefore, the 
flow variables in the leading-edge portion of the airfoil would have the 
values predicted by classical boundary-layer theory. This assumption 
removes the need to use Equation (16), the linking equation, in the 
airfoil stagnation point region where the thin-airfoil approximations 
are invalid. 

For ease of computation, the boundary-layer equations were 

14 

integrated using the Karman-Pohl hausen method as described in Schlichting . 
This integral method was chosen because it is accurate in regions of 
accelerated flow and also because the integral quantities it provides as 
output are just those needed in the interaction scheme. A more accurate 
finite-difference calculation could easily be used to replace the Karman- 
Pohl hausen scheme. 

The Karman-Pohl hausen method was applied to the Joukowski 

airfoil with the previously calculated u g distribution. The separation 

point was predicted at approximately s = .66S where S is the arclength of 

the airfoil. At tne stagnation point, the method yields 

* du 

a = 3.192, 5 = •641/— g^-| s _Q 
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It should be noted that the displacement thickness distribution 

agrees more closely with the finite-difference results than does the profile 

shape parameter a distribution. 

To begin the computational procedure, an initial estimate for 

tan o(s) is needed. Equation (15) shows that this requires initial 

* 

estimates of the 6 and <5 distributions. Since the tan 0 distribution will 

be continually updated in our interation procedure, a reasonable estimate 

is all that is necessary. The Karman-Pohl hausen results are used up to 

separation and then the estimates are guided by the experimental results 

15 

of Preston and Sweeting . For simplicity, the boundary layer thickness 

is assumed to increase linearly even though in the far wake it should become 

1/2 * 
proportional to s . tan 0 is much more sensitive to 6 than to 6 and some 

* 

numerical experimentation was necessary to obtain a reasonable 6 estimate. 

* 

6 was assumed to increase linearly from the separation point to the 
trailing edge and then to follow the equation 

6* = 1 + 12(S/s) 1/2 

in the wake which exhibits the proper asymptotic behavior in the far wake. 

A description of the computational procedure is given below. 

(1) tan 0 (s) is calculated from Equation (15) using Simpson's Rule for 

★ 

the numerical integration. <5 and 6 are held fixed at their boundary- 
layer values for the first 15% of the arclength so that tan 0 need only 
be calculated for s ^ .15S. The integration proceeds up to s = 15S 
since integration further downstream did not appear to affect the results 
on the airfoil . 
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(2) Equations (14) are integrated numerically in the interval 

.15S ^ s ^ 15S using a Runge-Kutta program from the University 

of Maryland Computer Science Center library, tan 0 is taken as the 

* 

average value at the endpoints of each interval. At s = .15S, a, 6 

and u g are taken from the Karman-Pohl hausen results. The integration 

* 

yields updated values for 6 ,6 , a and u g . 

•flf 

(3) The cycle is repeated with the newly calculated values of 6 and 6 
being used in (1) 

* 

Far downstream, u g should approach and 6 should approach 

a constant value. For a solution to be obtained the iteration procedure 

above must converge and be independent of the starting point for the 

strong interaction integration (taken arbitrarily above as .15S). 

It is difficult to assess whether the above procedure has in 

fact converged. It would appear to depend to some extent on which 

variables are being monitored. A typical computation is demonstrated 

in Table 1 where the variables shown are the separation point location, 

* * 

a, 6 and u g at s = 15S and <s at the airfoil trailing edge. The separation 

point moves closer to the airfoil leading edge with each iteration (eight 

iterations are recorded). The resulting solution is highly suspect since 

it seems unlikely that for the thin airfoil under consideration the 

separation point would occur so great a distance upstream of the value 

predicted by classical boundary-layer theory. 

* 

Some typical distributions of 5 , u g and tan 0 are shown in 
Figures (5-9) for three iteration cycles. Note that some of the variables 
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deviate immediately from the Karman-Pohl hausen results at the start of 
the strong interaction region. This is in violation of the assumption 
made regarding the limited upstream influence. Attempts to improve the 
convergence properties of the scheme by using under-relaxation in updating 
Equation ( 1 5) met with failure. 

Solutions were attempted with the starting point taken further 
from the leading edge and when the strong interaction computations were 
begun for s 2 .40S, the scheme appeared to converge but the computed 
separation point proved to then be dependent on the starting point. For 
example, for starting points of .4QS, .45S and .50S, the calculated 
separation points were .49S, .54S and .58S. This is clearly unacceptable. 

The results at the end of 10 iterations for the starting point of .50S 
are shown in Table 2. 

B. Weak Interaction Equations in Airfoil Forward Portion 

To be more consistent in the analysis, it is felt that the 
equations used in the forward portion of the airfoil should be the 
same as those used elsewhere. This should lead to smoother joining at 
the start of the strong interaction region. Consider the differential 
equations (9-11). In the leading-edge region, it is assumed that u g 
is given from the potential flow solution. Equation (9) is now eliminated 
and Equations (10-11) will be integrated to obtain 6* and a, using the 
known u e distribution. This abbreviated set is known as the weak interaction 
equations. 

* 1/2 

6, 6 and tan 0 are again scaled with v and the first derivatives 

are obtained in the form of Equations (14). 
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( 18 ) 


★ 

d6_ _ Njl da _ N? 
ds ” N * ds N 

where 

„1 - [(P " . R) ♦ 6* 2 (30 - (2H+1) ffi-)](gj/u e ) 


* * dU a 

N 2 = (RH-PJ)/6 U e + 6 -gf- J(l-H)/U e 


N = **3>air- J > 


The computational scheme now proceeds essentially the same as 

in Part A. The weak interaction equations are integrated from the leading 

edge to a point at the start of the strong interaction region, a distance 

far upstream of separation. This solution provides the starting values 
* 

of j , a and u fi for the strong interaction computation. The weak interaction 

equations predict separation at approximately 50% of the arclength so the 

* 

initial estimates for s and <5 must be suitably modified. 

Care must be taken in starting the weak interaction solution 

at the leading edge. Note that at the leading-edge stagnation point, 

* 

both 3s’ and are P ro P° r tional to u e • Since these derivatives must 
both be finite at s = 0, the following relationships can be obtained at 
s = 0 from Equations (18). 

3PJ = R(2H+1 ) (19) 

* du 

6 2 = P/(2H+1 ) -£\ s=q 
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The first expression leads to a = 2.967. The equations have 
a saddle point singularity at the leading edge. For the first .1% of the 
arclength, a similar solution is assumed so that a = 2.967 = constant. 
Equation (10) is integrated analytically from s = 0 to s = .Q01S using 
average values of u fi and du g /ds in the interval. The result is 


6* 2 (.001S) _ 
<5 * 2 ( 0 ) 


du e 
ds 


du. 


[2 -gg- (0) + (-g^- (.001S) 


du e 
• ~d? 


(0))e 


2(2H+1 ) 
H 



du Q 

-af-) ( .oois) 

avg j 


du e 

avg 

where avg stands for the average of the end point values. Equations (18) 
are numerically integrated downstream from s = .001S with care being taken 
to use small intervals in the leading-edge region. The weak interaction 

it 

solutions for a and 6 are compared to the Karman-Pohl hausen results in 
Figure 10. 


The complete scheme was run with the strong interaction region 
starting at both 15% and 30% of the arclength. There is no appreciable 
difference from the results of Part A. 

Another suggestion (Klineberg^) is to match the momentum thickness 
at the trailing edge rather than the velocity on the dividing streamline. 

ic ★ 

Since H = e/6 and 6 is continuous, this is equivalent to forcing H to 
be continuous at the trailing edge. In region III, 


H = .24710 - .43642a - .04773a 2 - .19654a 3 + ,41918a 1 * 


At the trailing edge, H is known from its value in region II 
and a is determined from the above equation. Use of this matching condition 
had little effect on the results. 
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C. Velocity Dependent Interaction Equation 

It lias been suggested (K1 ineberg^ 6 ) that the reason for the lack 
of success of the scheme is the possibility that Equation (15) and the 
integral continuity equation are not independent equations. Therefore, a 
new form of the equation linking the viscous and inviscid regions will 
be considered. It is obtained by a simple modification of the result given 
in Alber^ for the incompressible wake. 

Let 0 and V be the perturbation velocity components of the 
potential flow. Since they are harmonic conjugates, they are related by 
use of the Cauchy Integral Formula and the thin-airfoil approximation 
in the following equations. 


U(x.y) - I f 


V'(x.y) ■ - if 

This representation is formally equivalent to representing the 
effective body by a vortex distribution on the chordline. Now, 

U(x,y) = U + U(x,y) and U/U «1 


Therefore, 


<*" e hnvi$c1d * 


. i r (x-5)[l-U(5.0)/U.] ds 
" J-»(x-S) 2 + (f+v 1/2 i) 2 

On the surface, the limiting process of boundary-layer theory 
yields U( § ,0) = u e ( § ) . Using this and the previously developed relationship 
between (tan ©J-jnyjsc-jd an< * tan Q» tlie new 'f orm the linking equation 
becomes 
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(20) 


tan 


0 - *il 


(x-§ )[l-u e (§ )/Uj 
(x-§ )2+ (f-tv^ 2 S ) 2 


- • - *-’ /2 ^ 


■k 

Note that tan q now depends on u g instead of 5 . 

To obtain an initial estimate for the tan e distribution, 
several u g distributions were tried but it proved difficult to choose a 
u g distribution which would lead to a reasonable tan 0 . The next approach 
was to assume a tan 0 distribution directly for the first iteration cycle. 

It was possible to select tan 0 distributions which led to 
apparently reasonable values of u g and the other flow variables in 
the first iteration cycle but which yielded an unacceptable tan 0 distribution 
for the second cycle. To demonstrate the sensitivity of the scheme, 
consider the two tan 0 estimates. 



tan 0 = 5 sin* s/2S 

s < S 



= -2(S/s} 3/2 

s > S 

(21a) 

and 

tan 0 = 6 sin* s/2S 

s <. S 



= -2(S/s) 3/2 

s > S 

(21b) 


The u distributions from the first iteration cycle are shown 
e 

in Figure 11. The resulting tan 0 distributions are shown in Figure 12. 
In Figure 13, the tan 0 integral is separated into its airfoil and wake 
components. The term due to the airfoil slope is also shown. It can 
be seen that the major differences in the two cases stem from the wake 
contribution to the integral and that tan 0 is very sensitive to the wake 
results. Note also that the airfoil slope term is dominant over a large 
portion of the airfoil. It is difficult to obtain a tan 0 distribution 
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in the second cycle which is not negative for a substantial portion 
of the leading-edge region and this leads to unreasonable predictions 
from the viscous computation. With the above formulation, it was not 
possible to obtain reasonable results for more than two or three iteration 
cycles. 

In the previous formulation of the linking equation, there 

is no contribution to the integral from the flowfield upstream of the 

leading edge. Equation (20), however, has a contribution from this 

region since u g f in the upstream neighborhood of the leading edge. 

u fi was calculated from potential flow theory for x< Q,y=0. This portion 

of the u distribution is shown in Figure 14. Inclusion of these 
e 

results improves the tan 0 distribution on the airfoil but has a harmful 
effect in the wake region, tan 0 with the upstream effect included is 
compared to the result with the latter initial estimate in Equation (21) 
in Figure 15. 

At this time it was decided that the form of the linking equation 
given in Equation (20) was too sensitive to the flow quantities to be use- 
ful and that the most productive path to follow was to return to the 
formulation of Equation (15) and work towards developing a means to 
correctly update the solution in the leading edge region with each iteration 
cycle. 

It is possible that the decision to abandon the formulation 

of the linking equation presented here was a mistake. In their success- 

8 

ful treatment of a similar problem, Klineberg and Steger also encountered 
difficulty in the initial iterations due to erroneous wake results. 
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Their fix was to initially choose the downstream boundary of their 
computation region to be close to the airfoil trailing edge and to 
gradually move it downstream as the calculations proceeded. Also, 
they used rather strong under-relaxation to improve the scheme's stability. 
D. Updating of Solution in Airfoil Forward Portion 

A major deficiency of the method as it has been used is that 
the weak interaction solution is fixed for the forward portion of the 
airfoil up to the joining point with the strong interaction region. 

The flowfield near the leading edge is not allowed to adjust to the 
downstream changes in the displacement thickness distribution. This 
enforces a certain dependence of the final solution on the choice of a 
joining point - it is not clear, however, whether this alone has prevented 
a solution from being obtained. 

It is necessary to develop a method to calculate the potential 

flow surface speed for the flow past body plus displacement thickness 

which is more accurate than first-order thin-airfoil theory in the important 

leading-edge region. This method must be computationally easy to use 

or the purpose of originally using the thin-airfoil approximations will 

1 7 

be defeated. Von Karman had developed such a method for the solution 
for the potential flow past axisymmetric bodies. The method has been 
modified here for use in the two-dimensional case. 

The airfoil chordline is divided into constant length intervals 
and a source distribution is introduced with the source strength being 



constant for each interval. It is therefore necessary to calculate 
a finite number of source strengths. A major cause of error in the 
thin-airfoil approach is caused by the linearization and the transfer 
of the body boundary condition to the chordline. Both of these 
deficiencies are overcome fn Von Karman's method since the exact 
boundary condition is used and it is satisfied on the actual surface. 

The boundary condition is satisfied at the surface points corresponding 
to the midponts of the constant source intervals. The result is a 
set of linear equations for the unknown source strengths which is solved 
numerically. 

The details of the method are described in Appendix B. A 
comparison of the results for the Joukowski airfoil with the exact 
solution and the first-order thin-airfoil result for surface speed is 
given in Figure 16. Agreement with the exact solution is very good for 
distances greater than about .03S from the leading edge. The smallest 
constant source interval that could be used without introducing fluctuations 
into the surface speed was .01c where c is the chord length. 

This potential flow method is now incorporated into the computational 
scheme. To Insure continuity from Iteration to iteration, the results 
for the airfoil itself were used in the first iteration cycle. The 
scheme proceeds as before except that at the start of .each cycle, the 
potential flow past the body plus its displacement thickness from the 
previous cycle is computed and the weak interaction equations are integrated 
to provide the joining values at the start of the strong interaction 
region. The effective body is terminated at the now open trailing edge - 
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this is consistent with current techniques for "exactly" calculating 

the flow past two-dimensional bodies with displacement thickness 

1 8 

added (see, for example, Stevens, Goradia and Braden , where the 
effective airfoil is terminated closely behind the trailing edge). When 
the displacement thickness is added, the smallest source interval 
that can be used without smoothing the results is approximately .02c. 

The results of three iteration cycles are shown in Figures 
17-20. The effective body after one cycle is shown in Figure 17. 

The displacement thickness distribution is given in Figure 18 and u g 
is given in Figures 19 and 20, the former showing the very small change 
in the surface speed in the leading-edge region from the first to the second 
iteration. The addition of the updating of the leading-edge portion 
of the flowfield does not appear to appreciably affect the results. 

The trend is indicated in Table 3. An acceptable solution still cannot 
be obtained. 

As a final note, another deficiency of the method is discussed. 

For the strong interaction region, tan © is calculated from the linking 

integral equation and used as input in the integral continuity 

equation. This continuity equation is not satisfied in the weak interaction 

region. The question arises, in consideration of the smoothness of the 

joining, as to the value of tan © that is computed at the joining point 

by use of the continuity equation. For the previous set of computations, 

tan © was calculated at the joining point from both equations and the comparison 

appears in Table 4;. It is noted that the discrepancy between the values 

is less than 10%. 
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CHAPTER IV 


SOLUTION OF KLINEBERG AND STEGER 

Klineberg and Steger 8 have obtained solutions for subsonic and 
transonic separated flow past circular arc airfoils. The main difference 
between their computational technique and the one reported here is that 
they obtain a solution for the complete inviscid flowfield by use of a 
finite-difference relaxation scheme. Their method will now be described 
in the context of the previously reported procedures. 

A. Viscous Flow 

The system of equations which describes the flow in the viscous 

layer is the compressible flow equivalent to Equations (9-11). The four 

* 

unknowns are chosen to be u , v , 6 and H where v„ is used instead of 

e e e 

tan 0 and H is chosen as the profile shape parameter. 

The flow over the airfoil is divided into weak and strong 

interaction regions. In the weak interaction region, u g is obtained from 

the inviscid flow calculation and the equivalent of Equations (10-11) 

are solved by Runge-Kutta integration. The equivalent of Equation (9) 

is then integrated to obtain v g . A similarity solution is used to start 

the calculation at the leading edge. In the strong interaction region, 

* 

v e is given from the inviscid calculation and 6 , H and u g at the joining 

* 

point are given from the weak interaction solution. 6 , u g and H are 
taken to be continuous at the boundaries of the four distinct profile 
regions. > 

B. Inviscid Flow 

The flowfield outside the viscous region is considered to 
be irrotational and the transonic small disturbance approximation and 
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thin-airfoil boundary conditions are assumed. The governing partial 

differential equations are the continuity equation and the statement 

of the irrotationality of the flow in terms of the velocity components. 

Mixed boundary conditions are specified since the viscous computation 

yields v g in the weak interaction region and u g in the strong interaction 

region. Figure 21 shows a map of the flow regions. The equations are 

solved using the finite-difference relaxation scheme described in Steger 

g 

and Klineberg . 

C. Complete Interaction 

A complete interaction is calculated by alternately iterating 
the solutions to the viscous and inviscid equations. The sequence of 
computations is described as follows: 

(1) Solve the inviscid equations for the given airfoil. Compute the 
surface pressure distribution. 

(2) Solve the weak interaction equations from the leading edge to the 
joining point. For the biconvex airfoils used, the joining point 
is chosen upstream of the maximum thickness location at 40% of 
the chord. 

(3) Integrate the continuity equation to obtain v e . Downstream of the 
joining point, v g is arbitrarily specified in the first iteration 
cycle. 

(4) Integrate the strong interaction equations using the assumed v £ . 
Compute the u fi distribution. 

(5) Solve the inviscid equations with the mixed boundary conditions 
corresponding to v g in the weak interaction region and u 0 in the 
strong interaction region. Compute u g in the weak interaction 
region and v fi in the strong interaction region. 
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(6) Alternate between viscous and inviscid solutions until convergence 

is achieved. The boundary conditions for the current viscous/inviscid 
computation are provided by the previous inviscid/viscous computation. 

It is noted that the viscous and inviscid problems are linked 
through the surface speed in the strong interaction region. This technique 
is therefore related to the attempt to use Equation (20) as the linking 
equation in the previously reported research. It is recalled that this 
attempt was abandoned due to the sensitivity of the scheme to the flow 
variables in the initial stages of the computation. 

It is perhaps significant that Klineberg and Steger encountered 
similar difficulties. Until their scheme converges to a solution, the 
viscous layer integration diverges in the wake and generates unacceptable 
values of u g for use as boundary conditions for the inviscid problem. 

This difficulty is resolved in two ways. First, the location of the 
downstream boundary of the computation region is chosen close to the 
trailing edge initially and is moved downstream as the scheme progresses. 
Second, under-relaxation is employed in updating the inviscid flow boundary 
conditions in the wake. The above solution is aided by the use of an 
interactive computer graphics system. The discrepancy in tan 0 at the 
joining point that was inherent in the previous research is absent from 
the Klineberg and Steger technique since the continuity equation supplies 
tan 0 for the inviscid flow computation in the weak interaction region. 

Results were obtained for a range of Reynolds numbers and Mach 
numbers for 6 % and 12% thick biconvex airfoils. The results are independent 
of the initially chosen v e distribution and joining point location. For 
some cases, an approximation to simulate a turbulent wake was included. 
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A typical solution is shown in Figure 22 to illustrate the distributions 
★ 

of u , 5 and H. M is the free stream Mach number, Re„ is the chord 
Reynolds number, SEP indicates the separation point and RSP indicates the 
wake rear stagnation point. The final downstream boundary is located 
only 2.5 chordlengths downstream of the trailing edge. It is noted 
that previously the calculated distributions did not appear to approach 
asymptotic downstream values until much further downstream. 

Comparison with the experimental results of Collins (unpublished) 
is shown in Figure 23. Agreement is good when the turbulent wake 
approximation is used. The effects of Reynolds number and Mach number 
are demonstrated in Figures 24 and 25. In these figures, C p is the 
pressure coefficient. 
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CHAPTER V 


CURRENT RESEARCH IN VISCOUS-INVISCID INTERACTIONS 

At present, much work is being done in the area of viscous- 
inviscid interactions in flowfields with separation present. 

Research teams at various government, university and industrial 
laboratories are studying problems where the free stream is subsonic, 
transonic, supersonic and hypersonic, laminar and turbulent, and where 
the governing equations considered are either some form of the boundary- 
layer equations or the complete Navier-Stokes equations. In this chapter, 
an attempt is made to survey the current state-of-the-art and to 
essentially present a preview of what can be expected to appear in the 
literature in the next few years. This survey is by no means all 
inclusive. 

The only solution to a subsonic interaction problem including 
separation to appear in the literature apparently is the solution of 

O 

Klineberg and Steger of NASA-Ames Research Center which is discussed 
in detail in Chapter IV. Their method employs an integral form 
of the boundary-layer equations and uses a finite-difference relaxation 
technique to solve for the inviscid flowfield. Due to the computational 
difficulties encountered in using the integral approach in the viscous 
region, no attempt has been made to extend the results. Klineberg is 
currently studying under what conditions regular solutions to the 
boundary-layer equations which admit regions of separated flow exist. 

Armed with this information, the group at Ames may again attack sub- 
sonic and transonic interaction problems using finite-difference techniques 
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in the viscous as well as the in viscid region. 

A finite-difference technique to integrate the boundary-layer equations 
through regions of reverse flow is presented in Erdos, Baronti and 
Elzweig^ in connection with a study of transonic viscous flow with 
interaction. The authors question whether a solution of the direct 
problem (airfoil geometry given) by an iteration process between 
the solutions in the viscous and inviscid regions will converge and 
therefore treat the indirect problem. The airfoil contour is specified 
over the forward portion of the chord and the pressure distribution 
over the rearward portion. The inviscid flowfield is obtained by a 
finite-difference relaxation of the transonic small disturbance equation. 

The viscous equations are integrated using a Crank-Nicholson type 
implicit finite difference technique. In regions containing reverse 
flow, the parabolic equations take on an elliptic-like character since 
the solution at a point depends on downstream information. This is 
accounted for by use of upwind differencing in the reverse flow region. 

The viscous flow solution technique is verified by a solution for 
the incompressible flow past an elliptic cylinder using the measured 
pressure distribution as input. The theoretical results agree well 
with those from the experiment. The inviscid and viscous region solutions 
are not obtained simultaneously in this paper to generate a solution 
to the transonic interaction problem. 

The method presented in the above paper has not been extended by 
the authors who are at Advanced Technology Laboratories. Instead, interest 
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has been generated in the more practical turbulent transonic interaction 
problem. Baronti is currently working in this area and has developed 
an integral technique for the integration of the equations in the 
turbulent viscous region. An integral form of the turbulent kinetic 
energy equation is used to take account of the history of the boundary 
layer. The energy equation and the mean equation of motion describe 
the flow in the viscous region. The inviscid flowfield is handled by 
a finite-difference scheme which integrates the transonic small disturbance 
equation. The viscous region and inviscid region equations have not 
yet been integrated simultaneously so that a solution is not now available. 
Baronti is presently considering the solution of the indirect problem. 

Another group which is studying transonic interactions is led by 

R. Melnik of Grumman. The second-order accurate, implicit finite- 

19 

difference scheme of Keller is used to integrate the viscous equations. 
The boundary-layer equations are written as a set of first-order partial 
differential equations. Central differences are used in the streamwise 
direction to yield a block tri-diagonal system. To eliminate the pressure, 
the normal derivative of the streamwise momentum equation is used. This 
increases the order of the system and requires an extra boundary condition. 
This condition is obtained by a combination of the pressure-displacement 
thickness relationship at the viscous-layer edge and the wall compatibility 
condition which relates the pressure gradient to the normal wall shear 
stress gradient. Extrapolation is used in applying this condition. 

The method is used to study transonic free interactions. An 
asymptotic linearized solution is introduced at upstream infinity and 
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the system of equations is integrated downstream. The integration 

marches into the reverse flow region before numerical instability 

occurs. To obtain a solution for larger downstream distances, one can 

use upwind differencing as suggested above by Erdos et al^ or one can 

2U 

proceed as in Reyhner and Flugge-Lotz and eliminate the convective 
acceleration term in the reverse flow region. Currently, Melnik 
is developing an asymptotic solution valid in the flowfield far 
downstream in the reverse flow region to determine the proper downstream 
condition which the unique solution must approach. 

Melnik is also currently studying the subsonic (incompressible) 

interaction that occurs in the neighborhood of the trailing edge of a 

finite flat plate aligned parallel to a uniform stream. The above 

finite-difference technique is being used to integrate the boundary- 

21 

layer type equations developed by Stewartson for the inner viscous 
layer. The interaction proceeds as follows. First, the displacement 
thickness is prescribed and the boundary-layer equations are integrated 
and the pressure determined. The interaction equation relates the 
pressure to the displacement thickness in terms of a Cauchy integral 
and it is inverted to obtain an updated value of the displacement 
thickness. The cycle is then repeated. Upon the successful completion 
of this study, the angle of attack case will be considered. 

Also, Melnik is attempting to extend the approach of Stewartson 
to the turbulent interaction problem. He feels that two-dimensional 
and axi symmetric laminar and turbulent interactions can be handled by 
use of finite-difference techniques coupled with a careful analysis 
of the asymptotic behavior of the solution. 

M. Werle and co-workers at the University of Cincinnati and the 



Aerospace Research Laboratories have been studying supersonic and 

hypersonic interactions. The method for solving the laminar, two- 

dimensional supersonic interaction problem is described in Werle, 

22 

Polak and Bertke . An implicit finite-difference scheme is used 

in the viscous flow region. The governing equations are the boundary- 

layer equations with interaction. An analytical justification for use 

23 

of this model is given by Stewartson and Williams who performed an 
asymptotic analysis of the compressible Navier-Stokes equations for 

supersonic flow over a flat plate. The numerical technique is similar 

20 

to that of Reyhner and Flugge-Lotz but the equations are written in 
Levy-Lees type variables. In the region of reverse flow, at points 
with negative streamwise velocity the convective acceleration term which 
induces numerical instability is set equal to zero. If a major portion 
of the viscous layer has reverse flow, it is necessary to introduce 
artificially positive convective terms. A key to the success of the 
technique is the control of the continuity equation and its coupling 
to the momentum equation. The continuity equation is integrated from 
the viscous-layer edge to the wall and iteration is necessary to 
correct the assumed value of the edge normal velocity component. 

The interaction model makes use of linear theory to relate the 
pressure to the inviscid streamline deflection angle for the flow past 
the body plus its displacement thickness. The system is solved as 
follows. The problem is an initial value one and a "shooting" approach 
is used to correctly satisfy the downstream boundary condition. Self- 
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similar solutions are used at the upstream station. The energy equation 
is integrated using assumed profiles. The momentum and continuity 
equations are then integrated. The interaction equation is used to 
update the edge conditions and the cycle is repeated. Recently, Werle 
has developed a relaxation scheme to satisfy the downstream boundary 
condition which significantly reduces computation time. 

Results for the free interaction problem agree well with those 

23 

of Stewartson and Williams . Reasonable results are obtained for 
flows with shock-wave and ramp-induced separation bubbles when the 
inverse problem is considered. 

24 

The hypersonic interaction problem is treated in Dwoyer . Werle, 

25 

Vatsa and Bertke have taken a step towards solution of the three- 
dimensional interaction problem by solving for the flow past a swept 
compression ramp in a Mach 3 stream. The three-dimensional boundary- 
layer equations with constant cross flow are solved in Levy-Lees 
type variables. Linearized theory is used to relate boundary-layer growth 
to displacement surface angle. The solution technique follows that of 
the two-dimensional case. 

L. Olson and R. MacCormack of NASA-Ames Research Center are currently 
studying the solution of laminar subsonic interactions by the numerical 
integration of the Navier-Stokes equations. The method used is the time- 
dependent finite -difference technique of MacCormack . To test the 
method for subsonic viscous flows, solutions have been obtained for flow 
into an inlet and flow past a flat plate cascade with Reynolds number 
based on channel width varying from approximately 20 to 150. A paper 
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is in preparation in which comparisons are made with existing 
incompressible steady state solutions. 

At present, they are studying the subsonic interaction in the trail ing- 
edge region of a finite flat plate aligned parallel to the stream. The 
Reynolds number is of the order of 100 and the Mach number is less than 
.4. Upon the successful completion of this endeavor, many avenues of 
extension are open. Apparently, the next step will be to consider 
adding a turbulence model to the flat-plate interaction. It is felt 
that at present the cost of exploring high Reynolds number solutions 
is prohibitive. 



CHAPTER VI 


DISCUSSION AND CONCLUSIONS 

This report describes a method proposed to treat the problem of 
steady, two-dimensional, laminar, incompressible, high Reynolds 
number separated flow past thin airfoils. An integral description of 
the flowfield in the viscous region is used and the interaction between 
the inviscid and viscous flowfields is provided for by use of a thin- 
airfoil integral. In the light of the failure of this method to 
generate a solution, it is useful to compare it with the methods 
discussed in Chapters IV and V to help explain the lack of success and 
to suggest the path that future research might follow. 

Solutions to this or related problems with large Reynolds number 

by means of numerical solution of the Navier-Stokes equations are not 

available at present and do not appear to be forthcoming in the near 

future. This is unfortunate since these "exact" solutions could be 

used to test the validity of the approximate model equations which 

have been proposed. The methods discussed here use the boundary-layer 

equations with interaction to describe the flow in the viscous region. 

Some justification of this model for a supersonic mainstream and self- 

23 

induced separation is given in Stewartson and Williams but the asymptotic 
limit as Reynolds number goes to infinity is considered. The subsonic 
interaction is on even less of a firm footing since the most extensive 
research has been on the interaction at the trailing edge of a flat plate, 
a case without separation or reverse flow. It is clear that the need 
exists to provide further justification for the use of boundary-layer 
equations to study separated flow. 
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In the viscous flow region, either an integral or a finite- 
difference representation is used. In the pioneering work of Lees 
and co-workers, use of the integral form was motivated by the 
complexity of the computation and the absence of a suitable finite- 
difference scheme. According to Melnik and Werle, schemes to integrate 
into reverse flow regions are well in hand and it no longer appears 
to be necessary to use the integral approach. 

For the inviscid flow computation, two alternatives also exist. 

Most researchers consider a linearized interaction where the flow 
deflection angle or pressure is related to the slope of the displacement 
thickness. Klineberg and Steger and Baronti use a finite-difference 
integration for the inviscid region - this provides for a conceptually 
more direct matching at the viscous -layer edge. 

The method described here is somewhat unique. In the viscous 

g 

region, the approach of Klineberg and Steger is followed. However, 
in the inviscid region, the approach is similar to the integral 
representation of Melnik. This may be where the difficulty arises. 

The approach is the only one to use both the integral continuity equation 
and the flow angle-displacement thickness slope integral simultaneously. 

The equations have the same form and may not provide independent information. 

The major difficulty encountered in the use of an integral form 
of the linking equation, either the description given in I IB or the one 
in IIIC, for the subsonic interaction case with separation is the need 
to update the inviscid flowfield in the forward portion of the airfoil. 

This requires an auxiliary inviscid computation and negates to some 
extent the simplicity achieved by using an integral approach in the first 
place. For this reason, it appears that at the present time the most 
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effective way to handle the solution of the in viscid region 
equations is with a finite-difference relaxation scheme. It would 
be useful in treating a more general body to eliminate the "small 
disturbance" limitations of some of the current schemes. The 
essential features of the flow in the viscous region are captured 
by the use of an integral approach however increased accuracy of the 
solution and an improved knowledge of the details of the flowfield 
can be obtained by the use of a finite-difference technique in this 
regi on . 

In summary, a technique is not yet available for the 
solution of subsonic viscous-inviscid interactions with separation 
although much research is currently in progress. It is recommended 
that finite-difference schemes be used to calculate the flowfields 
in both the viscous and inviscid regions. 
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TABLE I 



A Typical 

Computation for 

the Method 

of IIIA 

Iteration 

Separation 


At s = 

15S * 

★ 

6 at 

Number 

Point 

a 

u e 

6 

Trailing edge 

1 

.47S 

.685 

.839 

6.93 

16.0 

2 

.38S 

.706 

.860 

6.28 

17.4 

3 

.34S 

.715 

.877 

5.97 

18.8 

4 

.31S 

.718 

.890 

5.82 

20.2 

5 

.30S 

.719 

.898 

5.75 

21.5 

6 

.29S 

.719 

.906 

5.69 

22.7 

7 

.28S 

.719 

.913 

5.67 

23.8 

8 

.27S 

.718 

.921 

5.67 

24.7 


Note: The starting point is s = .15S 

TABLE 2 

Computation for the Method of IIIA with Starting Point at .50S 


* 


Iteration 

Separation 


At s = 

• 15S * 

6 at 

Number 

Point 

a 

u e 

5 

Trailing edge 

1 

.62S 

.722 

.856 

5.96 

14.5 

2 

.60S 

.768 

.898 

4.73 

14.4 

3 

.59S 

.795 

.926 

4.07 

14.4 

4 

.58S 

.810 

.946 

3.71 

14.4 

5 

.58S 

.819 

.958 

3.49 

14.4 

6 

.58S 

.825 

.966 

3.37 

14.4 

7 

.58S 

.829 

.971 

3.28 

14.4 

8 

.58S 

.832 

.974 

3.22 

14.4 

9 

• 58S 

.834 

.973 

3.16 

14.4 

10 

.58S 

.837 

.975 

3.11 

14.4 



TABLE 3 


A Typical Computation for the Method of HID 


Iteration 

Separation 


At s = 

15S* 

6 at 

Number 

Point 

a 

u e 

6 

Trailing edge 

1 

.42S 

.668 

.850 

5.16 

10.7 

2 

.38S 

.703 

.883 

4.42 

11.2 

3 

.36S 

.723 

.903 

4.03 

11.7 


Note: The starting point is s = .15S 


TABLE 4 

A Comparison of tan 0 at the Joining Point Calculated from the Continuity 
and Linking Equations 


Iteration 

Number 


tan 0 from Equation (9), tan 0 from Equation (15), 
Integral Continuity Equation Linking Equation 


1 

2 

3 


2.043 

2.344 

2.492 


2.019 

2.125 

2.297 


Note: The computation is the same one used in Table 3. 
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APPENDIX A: POLYNOMIAL CURVE FITS OF PROFILE QUANTITIES 


Table A1 

7 k 

Coefficients of. Functions F =1 c k a for Regions I, II, III 

k=0 

F C Q C ] C 2 C 3 C, C 5 C g C ? 


Region I 


H 

0.24711 

0.11056 

J 

0.37372 

0.16969 

Z 

1 .03539 

0.48373 

R 

1 .25782 

-0.55550 

P 


0.48745 

dH/dA 

0.11056 

-0.04245 

dJ/dH 

1.50031 

0.28105 


H 

0.24711 

-0.25057 

J 

0.37372 

-0.42859 

Z 

1 .03539 

-1 .02605 

R 

1.25782 

1.09088 

. P 


-1 .19450 

dH/dA 

-0.25057 

-0.86024 

dJ/dH 

1.50031 

-0.84045 


H 

0.24710 -0.43642 

J 

0.37368 -0.59326 

Z 

1 .03285 

-1.35026 

R 

1.25775 

2.02922 


-0.02122 

0.00435 

-0.00097 

-0.02336 

0.00572 

-0.00175 

-0.01502 

0.02610 

-0.00370 

0.31964 

-0.09077 

0.01398 

-0.09927 

0.00960 

-0.00031 

0.01304 

-0.00389 

0.00050 

-0.04287 

0.00262 



Region II 


-0.43012 

0.1430 

-0.4267 

0.33036 

-5.1517 

10.5964 

-1.12405 

-1.1456 

3.3434 

7.01736 

-33.8762 

196.7688 

-0.70990 

-7.1253 

20.8568 

0.42888 

-1.7068 

-54.2937 

3.32376 

-13.8668 

5.4767 

Region III 


-0.04773 

-0.19654 

0.41918 

-0.00220 

0.27426 

-0.39421 

0.05127 

0.84914 

-1.89847 

3.88529 

-9.20873 

16.34366 


0.000099 

0.000191 

“0.000935 


-10.8587 38.7425 -31.209 

- 5.8174 

371.9762 244.3095 

100.2729 310.2394 -263.587 

232.4553 -218.4664 
30.1770 


dH/dA and dJ/da are obtained by differentiating the expressions for H and J, 
respectively. Then, dJ/dH * dJ/da/dH/da 



In Region IV, the following curve fits apply. 

6 = 2.35860 - 5.27379a + 6.73235a 2 - 6.12945a 3 + 2.34196a 4 

H = .58530 + .59289a - 2.14390a 2 + 1.48259a 3 - .52313a 4 

J = .88674 + .67519a - 1.66397a 2 + .28925a 3 - .19051a 4 

Z = 2.43208 + ,98139a - .91905a 2 + 1.01789a 3 - .52020a 4 

R = .53431 - .58670a - .62596a 2 + .89792a 3 - .22179a 4 

The profile quantities are 

H = H/?, J = I/s', Z = Z/S, R = RS 

Also, 

dH = £ dH/da - H d6/da 

3a " 5 2 

dJ SdJ/da - J d6/da 

3? p 

and 

gjj-= dJ/da/dH/da 
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APPENDIX B: VON KARMAN POTENTIAL FLOW METHOD 


Consider the potential flow of a uniform stream of speed 
past a two-dimensional symmetric body given by the equation y = f(x) . 

The disturbance caused by the body will be represented by N source elements 
placed along the chordline (x-axis), each with constant strength over an 
interval of length ax. The source with strength q.. per unit length lies 
in the interval extending from x. to x^. The velocity components of the 
flow induced by N of these sources at the point (x,y) are 


; q< Y 2 + (x-x.)2 

u (x.y> ■ l log - — - 

, y 2 + (x-x j+] 


V(x.y) * I ^ (tan- 1 ^ 


- tan 


1 


i+1 


1 _x 

x-x. 


The exact surface boundary condition is given by 


V(x,f(x)) 

U +U(x,f(x)) 


df 

dx 


(Bl) 


(B2) 


If the boundary condition is satisfied at N surface points 
( x k ,y k = substitution of Equation (Bl) into Equation (B2) leads, 

after some manipulation, to the equations 

j, A ik«i ■ - & < x k> k = 1 ' 2 - 3 N < B3 > 


where 


and 


Q. = q i /2U 00 w 

-1 y k 


A ik = tan " tan 


-1 y k 


+ \ gj’. ( x k ) log 


V x i + 1 


y k + ( J 


k- x i) 


y£ + (xk-x. 


k i+l) 
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Equation (B3) is a set of N linear equations for the N unknown 
values of the normalized source strength Q i . The equations can be put in 
matrix form and solved by a variety of techniques. Here they were solved by 
the Gauss-Jordan elimination method program from the University of Maryland 
Computer Science Center Library. 


Once the Q^'s are obtained, the surface speed on the body can 
be computed. The components of velocity on the surface, U k and V k , are 
determined with the help of Equation (Bl) and after some manipulation become 


■ u. *1 


¥ U »V_ y k + (x k - X 1 ): 


yl + (*i, - x.,,) : 


„ 8 „ + y2 k , 

V, — > y q. qoS L 1 — T7T ’ n 1 ' ■ ^ yd 

k i=i C< x k - x i> 2 +yg3 ,/2 [<y* i+1 ) 2+ y^ 1/2 


The surface speed is 

“e< x k> * < u2 k + V2 k> ,/2 


(B5) 


For the Joukowski airfoil, the equation of the surface is 
y k = e(2-x k )(2x k -x|) 1/2 

One would expect that the accuracy of the solution would increase 
with the increase in the number of constant source intervals. This is true; 
however, there is a limit on the number of intervals that can be used 
for a particular geometry. For a 10% Joukowski airfoil, the maximum number 
of intervals possible was approximately 100. When the number was increased 
beyond this, oscillations appeared in the singularity distribution and 
in the resulting surface speed. 

There is not complete freedom in the choice of the endpoints 
of the singularity distribution. For a closed airfoil, at the leading 
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and trailing edges .y^ = 0. To compute the surface speed at the leading 
edge, = 0 must be substituted into Equations (B4) and it is clear that 

if the singularity distribution begins at the airfoil leading edge and 
therefore x.. =0, becomes unbounded. For a particular geometry, there 
is a limitation in the distance between the leading edge and the start 
of the singularity distribution. For the 10% Joukowski airfoil, this 
distance was .5% of the chord. Starting the distribution any closer to the 
leading edge caused large oscillations in the surface speed in the leading- 
edge region. 

The boundary conditions are satisfied at the surface points 
corresponding to the midpoints of the source elements. The following 
placement of the surface points and the singularity points x^ is chosen. 
x k = kAx 
x i = ( i -1/2) Ax 

For the 10% Joukowski airfoil, x = .01c is chosen so that the first 
source element begins at x = .005c. 
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Figure 3. Coordinate System for Inviscid Region 



Figure 4. Joukowski Airfoil, Its Slope and Surface Speed 
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Figure 8. Typical Displacement Thickness Distribution for 
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Figure 10. Solution of Weak Interaction Equations for Profile 
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Figure 11. Inviscid Surface Speed Distributions Resulting 

from Initial tan 0 Distributions of Equation (21) 
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Figure 12. 



Updated tan 0 Distributions Resulting from Initial 


tan 0 Distributions of Equation (21) 
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Figure 14 . 


Potential Flow Surface Speed Distribution 
Upstream of Leading Edge 


63 



Figure 15. Updated tan e Distribution Resulting from 

Initial tan o Distribution of Equation (21b) 
Including Upstream Correction 
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Figure 17. Typical Airfoil Plus Displacement Thickness for 


Method of IIID 


Figure 18. Typical Displacement Thickness Distribution for 
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Figure 20. Typical Inviscid Surface Speed Distribution for 
Method of II1D (In Wake) 
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Figure 21. Computational Regions with Boundary Conditions 
for Method of Klineberg and Steger 
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Figure 22. Typical Distributions for Method of Klineberg 
and Steger 



Figure 23. Comparison of Pressure Coefficient with Experimental 
Values of Collins for Method of Klineberg and Steger 
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Figure 24. Effect of Reynolds dumber on Pressure Coefficient 
for Method of Klineberg and Steger 
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Figure 25. Effect of Mach Number on Pressure Coefficient for 
Method of Klineberg and Steger 
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